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The Landau-Zener transition and the surface hopping 
method for the 2D Dirac equation for graphene* 
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Abstract 

A Lagrangian surface hopping algorithm is implemented to study the two dimensional 
massless Dirac equation for Graphene with an electrostatic potential, in the semiclassical 
regime. In this problem, the crossing of the energy levels of the system at Dirac points 
requires a particular treatment in the algorithm in order to describe the quantum transition- 
characterized by the Landau-Zener probability- between different energy levels. We first derive 
the Landau-Zener probability for the underlying problem, then incorporate it into the surface 
hopping algorithm. We also show that different asymptotic models for this problem derived in 
|25| may give different transition probabilities. We conduct numerical experiments to compare 
the solutions to the Dirac equation, the surface hopping algorithm, and the asymptotic models 
of [25]. 

Keywords: Dirac equation; Wigner transform; Semiclassical model; Band crossing; Landau- 
Zener formula; Surface hopping algorithm; Spectral methods. 
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1 Introduction 

We are interested in the description of the transport of electrons in a single graphene layer. This 
material is a two-dimensional flat monolayer of carbon atoms which displays nnusual and interesting 
electronic properties arising from the bi-conically shaped Fermi surfaces near the Brillouin zone 
corners (Dirac points). The electrons propagate as massless Dirac Fermions moving with the Fermi 
velocity vp, which is 300 times smaller than the speed of light vp ~ 3 ^ ~ 10® and their 

behavior reproduces the physics of quantum electrodynamics but at much smaller energy scale. 
Although this model has been studied for a long time, see |5] for a bibliography, it has remained 
theoretical until the work of [27] where the graphene was produced for the hrst time. After these 
results, the interest of researchers on this material has shown a remarkable increase including 
applications in carbon-based electronic devices m and numerical simulations, see e.g. [^ and 
references therein. 

In this paper, we will consider a model of a two-dimensional Dirac equation ([T][3l[26]) of a 
graphene sheet in the presence of an external potential. This model consists of a small parameter 
h directly related to the Planck constant. We are interested is the design of an efficient numerical 
method-the surface hopping method- for the graphene Dirac equation in the semiclassical regime 
where h << I. In this regime, the solution of the Dirac equation is highly oscillatory thus a huge 
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computational cost is required to give accurate wave functions or physical observables for either 
finite difference methods ([a,[ii]) or time-splitting spectral method m, since one needs to resolve 
the high frequency both spatially and temporally. 

The development of efficient numerical methods for the related Schrodinger equation in the 
semiclassical regime has motivated many works in the last decade, see the review paper m and 
references therein. In the semiclassical regime, one often uses asymptotic analysis, such as the 
WKB analysis and the Wigner transform to help to reduce the computational costs and to develop 
efficient computational methods based on the asymptotic models. In the framework of Wigner 
transform, the idea is to construct a measure on the phase space, called the Wigner measure, when 
/i —>■ 0, to obtain the physical observables (such as density, flux, and energy) with a computational 
cost far less than a direct quantum simulation. When the gap between different energy levels is of 
order one (the so-called adiabatic case), the Wigner measure technique provides a simple description 
of the motion; it can be well approximated by a fully diagonalized system, one classical Liouville 
equation for each energy level m- However, in the graphene Dirac equation, the energy levels 
cross at the Dirac points, where non-adiabatic transfers are observed and the particles can tunnel 
from one band to the other. The standard Wigner approach then needs to be revised to describe 
the non-adiabatic phenomena. 

One of the widely used approaches to simulate the non-adiabatic dynamics is the surface hop¬ 
ping method initially proposed by Tully and Preston [34] as an efficient computational method to 
go beyond the classical Born-Oppenheimer approximation. This method is widely used in chem¬ 
istry and molecular dynamics, see for examples [3 [301 [331134]- The basic idea is to combine the 
classical transports of the system on individual potential energy surfaces with instantaneous tran¬ 
sitions from one energy surface to the other. The transition rates for this band-to-band hopping 
are given by the well known Landau-Zener formula |35| . From the mathematical point of view, the 
first rigorous analysis of non-adiabatic transfer using Landau-Zener formula dates back to Hage- 
dorn m- More recently, the Wigner measure techniques for separated energy levels have been 
extended in [7] and [8] to systems presenting band crossing. The proof is based on microlocal 
analysis to justify the Landau-Zener formula. Departing from the results in [7] and [5], a rigorous 
surface hopping algorithm in the Wigner picture was proposed in [5D1 for time-dependent two-level 
Schrodinger systems with conically intersecting eigenvalues and implemented numerically in the 
Lagrangian formulation in m- The corresponding Eulerian numerical scheme was proposed in 
|18) by formulating the hopping mechanism as an interface condition which is then built into the 
numerical flux for solving the underlying Liouville equation for each energy level. 

In the present article we give a Lagrangian surface hopping algorithm for the graphene Dirac 
equation similar to the algorithm in m- First, it is a classical result that the Wigner trans¬ 
form leads to two decoupled classical Liouville equations for each energy level, under the adiabatic 
assumption [lOj . At the Dirac points where non-adiabatic transition occurs, we first derive the 
Landau-Zener formula, and then incorporate it into the surface hopping algorithm. We also show 
that a reduced asymptotic model developed in |25| could give incorrect transition probability. We 
then compare through several numerical examples the solutions of the Dirac equation (solved by 
the time-splitting spectral method), the surface hopping algorithm and the asymptotic models of 
|25) . Our numerical results show that, when there is no wave interference, the surface hopping 
algorithm indeed gives the correct non-adiabatic transition at Dirac points, with a much greater 
computational efficiency compared with the simulations based on solving directly the Dirac equa¬ 
tion. 

The article is organized as follows: in section [2| we give the graphene Dirac equation and 
its semiclassical limit via the Wigner transform in the adiabatic case. We give some examples 
of potentials to show that non-adiabatic transition is indeed possible. In section [3 we derive 
the Landau-Zener transition probability for the graphene Dirac equation. The surface hopping 
algorithm is given in section 0] In section [3 we study the asymptotic models introduced in 
and show that a reduced model could give the incorrect transition probability. Numerical results 
are given in section [3 for comparisons of different models. For reader’s convenience we give the 
time-splitting spectral method of the Dirac equation in Appendix [^ 
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2 Quantum transport in graphene and the Wigner measure 


We consider the description of the transport of electrons in a single graphene layer in the presence 
of an external potential. Following [nisi US], it is modelled by the two dimensional Dirac equation: 

J ihdti} = [—ihvpo'idxi — ihvFO'2dx2 + 9^] V'j t G M, x G ... 

\ ip^Oyx) = -tpiix), X G 

where h is the reduced Planck constant, vp the Fermi velocity, q is the elementary charge of the 
electron and V(x) G M is the electric potential. The Pauli matrices cti , tT 2 are given by: 



The initial wave function ipiix) G is normalized such that: 



\tj}i{x)\'^dx 


1 


and, using the mass conservation, the wave function 'tjj{t,x) G satishes 


( 2 . 2 ) 


/ \tp{t,x)\'^dx = 1 (2.3) 

where t, x = (a:i, X 2 ) are the time and space variables respectively. We consider the system (12.11) in 
the semiclassical regime. For this purpose, we rewrite the equations such that there remains only 
one dimensionless parameter h. Proceeding as in ED, we change the variables to dimensionless 
variables as follows 

t^t/To, x^x/L (2-4) 

and define 

ui{x) = Li(ji(Lx), u{t,x) = L'ijj(Tot, Lx) (2-5) 

where T is a reference length and Tq = L/vp. We remark that ui(x) and u(t,x) are chosen such 

that the change of variable preserves the normalization (12.21) and (12.3|) . Plugging (12.51) into (EH), 
and dividing by the reference energy mvp where m is the effective mass of the electron, one gets 
the dimensionless graphene Dirac equation: 


J ihdtu = [—ihaidxi — iha 2 dx 2 +U]u, t G IR, x G 
\ u{0,x) = ui{x), X G IB?' 


where 

h = ^ 

mvpL 

is a small dimensionless parameter and U(x) is the dimensionless potential dehned by 


(2.7) 


U{x) = -^V(Lx). 

mvp 

We will consider the solution of (12.61) in the limit h —>■ 0. 

Non-adiabatic transfer happens at the Dirac points which are the crossing points of the eigenval¬ 
ues of the symbol related to (12.61) . To be more precise, we consider the complex 2 x 2-matrix-valued 
symbol: 

Po{x, 0 = B{0 + U(x)I, (x, JRl , (2.8) 

where I is the 2x2 identity matrix, B{^) is given by 
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and U S C°°{1R'^,1R) is such that V/3 € INq there is a constant > 0 verifying: 

\d^U{x)\<C^, VxGfR 2 _ 

An easy computation shows that the matrix B(^) has eigenvalues ±|^| with corresponding or¬ 
thonormal set of eigenvectors given by 

Moreover, it holds Pq G S^{1R'^ x and the Dirac equation (12.61) can be rewritten as: 

f ihdtu^ = Po{x,hD)u'', tGR 12 10) 

\ u^{0) = Uj 

where Po{x, hD) is the Weyl operator defined for u G by the integral: 

Po{x,hD)u{x) =-^ [ f Po u{y)P^^~y^-^d^dy, (2.11) 

(27r)^ 7 r 2 7jfj2 V 2 J 

and D = —idx- The operator Po{x, hD) is essentially self-adjoint on the Hilbert space TL = 

and the domain of its self-adjoint extension is Therefore —j^Po{x,hD) generates a 

strongly continuous group of unitary operators solution to ( 12 . 101 ) . 

The eigenvalues X±{x,^) = U{x) ± |^| of the symbol Po{x,^) satisfy A+(x,^) = X-{x,^) at the 
crossing set {^ = 0} C ^ The semiclassical limit away from the crossing set was performed 
in m for systems of the form ( 12 . 101 ) with more general symbols and initial data u’} subject to 
additional conditions. In m, the authors show that the semiclassical limit for the different bands 
can be treated separately where, for each band, the related eigenvalue plays the role of a scalar 
classical Hamiltonian. 


2.1 Adiabatic semiclassical limit 


We recall now some basic notions of the Wigner analysis involved in the semiclassical limit. For 
u, V G L‘^{1R'^), the Wigner transform is defined by: 

w^{u, v){x, ^) = J^^u(x- u -f P^^dy 

and for u G "H, the 2x2 Wigner matrix is defined by: 

We denote by w^[u] = tTW^[u] the scalar Wigner transform of u. For any bounded sequence 
in H, there is a subsequence of W^[f^] which converges in S'. Such a limit is called a Wigner 
measure associated to f^. If admits only one Wigner measure, we shall denote it by IF°[f^l 
and set w°[f^]=trW°[f^]. 

Another important object is the classical flow 4>^{x,^) = {x^(t), (t)) corresponding to the 

eigenvalues A±, i.e. the solution to: 


= x±(0)=x 

i^^{t) = -dxU{x^{t)), e^(0) = C 


( 2 . 12 ) 


where (x,^) G IR^ ^ Indeed, the decoupled semiclassical limit is valid on a set of the phase 
space which is stable under the flow (j)^ and where no band crossing occurs. More precisely, if 
there exists an open subset H C x IR| such that: 

H D {^ = 0} = 0 and cjif (H) C H, Vt G IR 
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and if the initial condition Uj in (12.1011 has a Wigner measure W'j such that uPj = tr satisfies 

w?|n<= = 0 

then using the results in [10] , it holds that [u ^], the scalar Wigner transform of the solution 
to ( 12 . 101 ) . converges to 

+w°_{t,x,i) (2.13) 

where .,.) is the scalar positive measure on x JR? solving the classical Liouville equations: 


dtw^ ± ill 9xW± - dxU.d^w^ = 0, IRt X n 




(o,.,.) = tr (n±wo), n 


= teM 


(2.14) 


For all {x,^) £ O, n±(^) in (12.141) denotes the projection on the eigenspace related to the eigenvalue 
^±{x,0- Using formula (12.91) . it can be expressed explicitly as follows: 


n±(e) = -(/±^i?(e) 


(2.15) 


Moreover, the density n^{t,x) = |u^(t,x)p converges to 


n°{t,x) = ( w°{t,x,^)d^. 

Jr'^ 

It follows that for initial data Uj such that the bands are separated initially, i.e. the support of 
does not intersect with the crossing set = 0 }, the measure is described by if the support 
of rcj is stable under the characteristics solution to (I2.12|) (in that case, the characteristics starting 
from the support of do not reach = 0 }). 


2.2 Some case studies 
2.2.1 Case U = 0 

In the case of the trivial potential, the solutions to (12.121) are given by: 

(a;=^(t),C=^(t)) = (2.16) 

for ^ 7 ^ 0 and we can take ft = IR^ x IR| \ {^ = 0}. Then the system (12.141) writes: 

I ± = 0 , Mt X 0 } ^ 

( w±(0,= Wj j_, {'f 7 ^ 0} ■u;^(t,x, 0) = 0, t £ M, X G 

where 

= tr (n±W°) , e^O. 

Using (12.161) in (I2.17L we get that: 

w±{t,x,^) = (^XT, Ct^O. 

Therefore, the density n^{t,x) = \u^{t,x)\'^ converges to 

„•«, X) = <_ (x + Xi,«) , (x - ijf. f) d{ . 

In the present case, if the support of the initial scalar Wigner measure does not contain the 
point {^ = 0 }, then no hopping occurs: the bands do not communicate and at any time the measure 
u>° is described by the separate evolution of the level characteristics. 
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2 . 2.2 Case U = axi, a G ]R \ {0} 

In that case, the solutions to (12.121) are such that: 


(2.18) 


^^(t) =^ - iat,0) 

and we can take = IR^ x \ {^2 = 0}. Then, system (12.141) becomes: 
f dtw^ ± = 0, Mt x {^2 ^ 0} , 

1 <(0,.,.) = tr (n±lT0), U 2 y^ 0 }, <(t,x, (^i,0)) = 0, t G M, x G G JR 

(2.19) 

Unlike in section 12.2.11 non-adiabatic transfer may occur. Indeed, if the support of does not 
intersect with the crossing set = 0 } but contains points of the form {x, ^) where ^ = (Ci, 0 ) and 
^ 0 , then the characteristic curve C^(t), given by (12.181) . starting from (^i, 0 ), will reach the 
crossing set at the time t = ^ and a band-to-band transition will take place. 

In these conditions, the system (12.191) does not describe correctly the asymptotics h —0 of the 
Wigner matrix W^[u^]. 

The goal of this paper is to develop efficient semiclassical methods to compute the non-adiabatic 
transition between different bands. To quantify the transition rates, we propose in section [3] a 
Landau-Zener formula which can be justified theoretically by using the results in [5]. Indeed, in 
[8] the following symbol was considered: 


Pl{x,0 

= A{0 + U{x)I, {x,0&PilxPih 

( 2 . 20 ) 

where / is the 2 x 2 identity matrix and 



- ( 1 -1 ) ^ 

( 2 . 21 ) 

Using the unitary equivalence 

BiO = R*A{OR, 

( 2 . 22 ) 

where 

* -* ) ’ 

(2.23) 


it follows that the two symbols have the same eigenvalues and Pq can be reduced to Pi after 
conjugation with the matrix R in the treatment of band crossing. 


Remark 2.1. The functions w^{t,x,^) in (12.131) are the diagonal terms of the Wigner measure 
W^{t,x,^) := lU°[u^(t)](x, ^). Indeed, it was shown in [TU] that is given by 

< = tr (n±lU°) In (2.24) 

and that is diagonal in the sense that = n+lU°n+ -|- n_lU°n_ on Mt x fl. Moreover, an 
easy computation leads to: 

n±w°n± = tr (n±iu°) n± = wln± on jr* x u. (2.25) 

3 The Landau-Zener transition 

In this section we will study the Landau-Zener transition for the Hamiltonian Po{x^ hD). 
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3.1 Classical flow around the crossing set 


As remarked in section [21 non-adiabatic transfer happens only when the characteristics reach the 
crossing set = 0}. Proposition 13.II below says that such characteristics will exist as soon as the 
potential has points x £ IR^ such that dxU{x) ^ 0. We refer to [5] for the proof. 


Proposition 3.1. Consider x £ IR^ such that dxU{x) ^ 0, then there exist two unique curves 
s I—>• (x^(s), ^^(s)) which are continuous in s in a neighborhood of 0 and for s ^ 0 and such 
that: 


j^x±(s) = ±i|^, x±( 0 )=x 

fs^^{s) = -dxU{x^{s)), e^( 0)=0 


To illustrate the fact that a spectral transfer for the symbol Pq happens at the crossing set, 
we will come back momentarily to the potential U = axi introduced in section 12.2.21 For such a 
potential, problem writes: 

I = a;±(0) = a; 

[ = J )> e^(o) = o ’ 

and \/x £ IR^, its unique solution is given by: 

x^{s) =XT sgn{a) ^ ( 0 ) ’ 


Plugging this solution in the projectors defined in (I2.15|) . we get for s < 0: 


n+(?±(s)) = - 


and for s > 0: 


n+(^±(s)) = - 


1 sgn(a) 
sgn(a) 1 


1 —sgn(Qf) 

—sgn(Q;) 1 


n_(e±(s)) = i 


-sgn(a) 


n-(e^(s)) = :r 


-sgn(a) 

1 


1 

sgn(Q!) 


sgn(a) 

1 


It is easy to see that the projectors n+(^^(s)) and n_(^=*^(s)) interchange one with the other when 
the characteristics pass through the crossing set = 0}. 


3.2 A heuristic derivation of the Landau-Zener formula 

In this section, we give a heuristic argument, similar to the one in m and m, to derive the 
Landau-Zener formula for the Hamiltonian Po{x, hD). 

In general, the region for non-adiabatic transfer is not restricted to the crossing set = 0}, 
since quantum transition occurs as long as the two energy levels are sufficiently close, in the case 
of avoided crossing [TT]. As in [T2] and [20] i we define this region to contain the points where the 
distance between the eigenvalues A±(a;,.f) ofPo(a:iC) is minimal. In our case |A+(a;,^) —A_(x,^)| = 
2 |^| and, when considered along the characteristics solution to (12.121) . the necessary condition for 
minimal gap is: 

(l$(s)P)' = 0 ^(s).a,,17(j:(s)) = 0 

and the hopping surface is chosen as the set: 

S={{x,O&lR^;^-dxUix)=0}. (3.2) 

The heuristics follows by inserting the characteristics in the trace-free part of the symbol to obtain 
the system of ordinary differential equations: 

ihfjfs) = P(^(s))V'(s). 
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After conjugation with the matrix R defined in (12.2311 and using equation (I2.22|) . one arrives at the 
following new system: 

ihij)'{s) = A(^(s))'0 (s). 

Assume the particles defined by the trajectory ( (12.121) 1 are near a point G S (due to 

translation invariance we assume particles are at initially), then the Taylor expansion 

gives: 

x{s) =x* ± ||^s + e>(s2) 

m=C-d.U{x*)s + Ois^). 

Ignoring the O(s^) terms, the system becomes: 

\ £.*2-dx2U{x*)s -Ci+d^^U{x*)s ) ^ 


After conjugation with the rotation matrix: 


cos 9 sin 9 
— sin 9 cos 9 


where 9 is the angle such that 


we get: 


(cos 20, sin 20) = 


d:,U{x*) 

\d.U{x*)\ 


\dxU{x*)\ 




—s 

/\dxU{x*) 


AdxUjx*) \ 

|9x(7(^*)P I 


|axC/(x *)|2 

where ^ A( = ^ 2(1 — C 1 C 2 for C G ■ After conjugation with the matrix 


0 1 
1 0 


it follows: 


^^xU{x*)\ ^ 

If we set £ = |g 0 ~ system becomes: 


- e/\d^u(x') 


i* f\dxU{,x*) \ 
|axC/(x -)|2 


ietl;'{s) = 


5 7] 

7] —S 


i’is) 


which is the well known Landau-Zener problem (see |35| 1 for which the transition probability is: 

T = e"?’'" . 

This allows us to propose the following non-adiabatic transition rate at the point (a;*,C)' 


„ («*ASxr/(x*))2 
■ h |0x,t/(x*)|3 


T{x\C) = e 


(3.3) 


















3.3 About the rigorous justification of the Landau-Zener formula 


As already noticed in section [21 the symbol Pq involved in the Dirac equation (j2.10l) satisfies the 
identity 

Piix.O = RPoix.OR* 

where Pi and R are given respectively by (I2.20|) and (12.231) . Therefore, if denotes the solution 
to (12.101) . the function = Ru^ satisfies the equation: 

ihdtv^ = Pi{x,hD)v^ (3.4) 

where Pi{x,hD) is defined as in (12.111) . A Landau-Zener formula was obtained rigorously in [5] 
for the two-scale Wigner measure iy of the function (see [5], [20] for the definition of two-scale 
Wigner measures). This result provides a rigorous proof of our Landau-Zener formula. Indeed, if 
lyj denotes the two-scale Wigner measure of u^, it follows that the two-scale Wigner measure of 
is 1 / = RvjR*. Then, the Landau-Zener formula for vi can be deduced from the Landau-Zener 
formula for v. 

Remark 3.2. The Landau-Zener formula obtained in |5] writes 

— 

T = e 


where 


rj = A 


d£_ 
\d^U\ ■ 


When the direction 5^ is equal to this corresponds to the Landau-Zener formula (13.31) that we 
obtained heuristically in section [221 


4 A surface hopping algorithm 

We give here a semiclassical Lagrangian algorithm for the evolution of the diagonal terms of the 
Wigner matrix W^{t,x,^) of the solution 'u^(t) to (12.101) . The algorithm is adopted from the 
method proposed in m for time-dependent two-level Schrodinger systems with conically inter¬ 
secting eigenvalues. 

Define the level populations: 

P^{t) = \\U^ihD)u^{t)\\l, (4.1) 

where for u £ R, Il±{hD)u is the function defined via its Fourier transform Il±{h^)u{^) and 
it = Tu is given by: 

Pu{^) = 7 ^ / u{x)e~^^'^dx 

(it is clear from (12.151) that Il±{hD)u = n±(P)tt). With similar computations as in [TO], we obtain: 

P±(.t)= [ [ w^{t,x,^)dxd^, (4.2) 

J]Rl j Rl 

where 

w\{t,x,C) =R (n±(C)W^^(^,a:,0) • (4.3) 

As it appears from the following equation 

n±VF'‘n± = , 

which is obtained in a similar way as (12.251) . are the diagonal terms of the Wigner matrix W^. 
Up to a small remainder, the function w!^ can be written in terms of the scalar Wigner transform 
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of the level projections of the solution u^{t). Indeed, using Lemma 2.3 in [TU], it holds for all 
t G M: 

w’l{t)=w>^[u>lit)]+oil) (4.4) 

in V X (IR| \ {^ = 0})^ when /i —>■ 0, where u^{t) is the function with Fourier transform 
given by 

«±(i,0 = n±(he)u^t,e)i^^o 

and u^{t) is the Fourier transform of u^{t). Comparing the relations (14.3|) and (I2.24L it follows 
that, in the situation of section [5] without band-crossing, the partial differential equation in (12.141) 
satisfied by can be used to approximate the time evolution of . In the case of band crossing, 
the idea is to use (12.141) for the time evolution of lu^ as long as the classical trajectories solution 
to (12.121) are away from the hopping surface S defined by (13.21) . When a trajectory reaches a 
point (x* ,^*) G S a non-adiabatic transfer of weight occurs between and with transition 
probability given by (13.31) . 


The algorithm 

1. Initial sampling: in this step, an appropriate sampling of the function Wj^ defined in (14.81) 
is chosen. Specifically, a set of sampling points 

{ixk,^kjk) G Ml X Ml X {-,+}; k = l,...,iV} 

are chosen with associated weights Wk G M given by: 

= ^Gjk^Xk.^k) ■ 


2. Hopping transport: away from the set S', each particle {xk,^k,jk) is transported by the 
associated classical flow solution to (12.121) . In other words, for t > 0 small enough: 

(a:fc(i),'ffc(i)) = (Af (a:fc(0),^fc(0)), Wfe(t) = Wfc(0). (4.5) 

If there exists t* > 0 such that {xk{t*),^k(t*)) =: {xl,^l) G S, the weight is reduced using 
the transition rate 

T* = T{xl,Ck) , 

where T{x,^) is given by (13.31) . Moreover, in order to describe completely the non-adiabatic 
transfer, a new particle with index I > N is created. Specifically, for t > t* 

{xk{t),^k{t)) = (a;fe(0),^fc(0)), Wk{t) = (1 - T*)wk{t*) 

and the new particle is created for t > t* 

{xi{t),^i{t)) = (l>f_t,{x*k,^l), ji = -jk, wi{t) = T*Wk{t*). 

3. Final reconstruction: at the final time t/ > 0, there are M > N points 

{{xk,^k,jk) G Ml X Ml X {-,+}; k = 1,...,M} 

with associated weights Wk which are approximations to {tf, Xk, Cfe)- Then, using equation 
(SJ), the surface hopping approximations ±{tf) of the level populations P±{tf) are given 
by: 

M 

T‘i)i,±(i/) = '^WkSf^uJk (4.6) 

fc=i 

where (5) is the Kronecker symbol related to i and j, and ujk is an appropriate quadrature 
weight. 
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Remark 4.1. We note that this surface hopping algorithm is subject to some restrictions. First 
only the dynamics of the diagonal components of the Wigner matrix away from the crossing set are 
well approximated. Second, there are possible interferences which are not captured if no particular 
treatment is performed. Indeed, if w^{t, x,^) and w^{t, x,^) arrive at the same time at some point 
{x, close to the crossing set then a transfer of weight using only the Landau-Zener formula (13.31) 
might give an incorrect approximation of the dynamics. To avoid these interferences, we make the 
assumption that the initial data Uj in ( 12 . 101 ) satisfies 

n_(hO*?(e) = o (4.7) 

where u’} is the Fourier transform of u^. Then, it follows from (14.41) that the diagonal terms of the 
initial Wigner matrix 

(x, 0 = tr (n± [4] (x, 0) (4.8) 

satisfy in V x (IR| \ {^ = 0 })^ 

Wj _^_ = [uj] + o(l), Wj _ = o(l) (4.9) 

when h —>■ 0. Therefore, for the potential U{x) = axi given in section [2.2.21 the condition (14.71) 
insures that no interferences occur. Indeed, if w'^{t,x,^) reaches the crossing set = 0} at some 
time t* > 0 , the particle corresponding to w^{t,x,^) immediately moves away as it appears from 
the equation 

^±(t) =-a(t-t*,0) (4.10) 

for the momentum part of the characteristics solution to (12.121) . Similarly, in the case where the 
potential has no stationary points, i.e. dxU{x) 7 ^ 0, Va: £ equation (14.101) is replaced by 

^^{t) = ~ (a;'^(s)) ds 

and the condition (I4T1) insures that no interferences occur as long as — is small enough. 

Remark 4.2. To deal with this interference, a possible solution might be to use a hybrid method 
for the Dirac eouation (12.101). Such a method was introduced in m for the Schrodinger equation 
and mixes a Gaussian beam method or a Liouville equation in the region where no-interferences 
occur and a complete quantum solver in the region where the phase information of the wave 
function is required. In our case, the second region corresponds to the hopping region. Since the 
region which involves the quantum solver can be chosen very small, a hybrid method allows an 
adapted treatment of the interferences without increasing dramatically the numerical cost even if 
the semiclassical parameter tends to 0. In addition, the resulting algorithm is supposed to work 
regardless of the conditions in Remark |4.II on the potential and the initial data. Another possible 
solution is to keep the off-diagonal entries, which contain information about the non-adiabatic 
transition in the Wigner transform, and then derive the semiclassical models for the entire Wigner 
matrix, see [25l|4]. In the next section, we study such a model. 


5 


Transition rate of an asymptotic model derived in 


25 


As it was mentioned, the surface hopping method breaks down where there are interferences. A 
possible remedy for this is to derive improved models using the Wigner transform for the en¬ 
tire Wigner matrix, in which the off-diagonal entries contain non-adiabatic transition information 
ESlll]. In this section, we study such a model obtained in m. and compare the Landau-Zener 
transition rate of this model with the one we derived in section [S] To serve our purpose we will 
only consider the model for the potential studied in 12 .2.21 U = ax\, ex £ IR \ {0}, for which the 
Wigner matrix does not depend on the variable xi. It was shown in |25j how this asymptotic model 
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can describe non-adiabatic transfer as a quantum correction of the decoupled system (12.1911 . 


When the Wigner matrix does not depend on the variable 3 : 2 , the system (12.191) reads 


f i9tw± ± = 0, IRj x {^2 0} 

1 <(0,.,.)=tr (n±W0), U 2 ?^ 0 }; x, (^i,0)) = 0, t G M, x G G M 

(5.1) 

By taking into account the scaling (12.41)(I2.5|l , the asymptotic model in is the following correction 

of dSH); 

I atw± ± - a%w± = ±a||^Im , Mt x {^2 ^ 0} _ 

\ dtWi - iA{^)wi - ad^^Wi = -if - w-) 


where 


m = 


h icp ■ 


In the limit h —>■ 0, the function Wi(t,x,^) approximates the off-diagonal terms of the Wigner 
matrix of the solution to (12.1011 . As it appears from the last equation in (EH), the function Wi 
depends on w+ and W-. As a consequence, wt provides a coupling term at the r.h.s. of the first 
two equations in (15.21) . 


Using the method of characteristics, Wi can be expressed as an explicit function of the difference 
w+ — W-. Inserting this solution in the first two equations in (15.2L the following approximate 
equations for ri;±, in the limit h —^ 0, are obtained in [25) : 

dtw± ± - ad(_^w± = =ft(^)(w+ - w_), Rt x Rx^ x (5.3) 

where ^2 0 can be considered as a small parameter. In the domain 

G R^ such that 0 < |^ 2 | < and |^i| < > 


the coefficient r(^) is given by: 


t{£.) = f (^|sgn(a6) - arctan^^ (5.4) 

where sgn(x) denotes the sign of x. As it will be the case in section 16.21 for our surface hopping 
algorithm, the set {^1 = 0} is the important region for the hopping. Indeed, the following lemma 
(proof left to the reader) shows that the set {^1 =0} plays the role of an interface where, in the 
limit ^2 —)> 0, the solution to (15.311 will have a discontinuity. 

Lemma 5.1. The function ^1 1 —> ^(^ 1 ,^ 2 ) tends in T)'[R) to |a|/?^jj=o when ^2 0, where 

TT^ 

/3=—. (5.5) 


In the case a > 0, we show below that non-adiabatic transfer is possible using the model (15.311 
and give the corresponding transmission matrix at the interface. Define uj± as the set: 

uj± = {Rt X Rxi X R^i) l~l {±^1 > 0} 

and consider an initial condition in the upper band and localized at the right of the interface. In 
other words, the initial conditions for (15.311 are 

■*^±(0,.,.) = 

where ic/,- = 0 and wi_^ is a function in C^{Rx^ x R^^) which is independent of ^2 and has 
support in uj+. To perform the limit ^2 0, the following assumptions are required on the solution 

to (15.31) : 


12 










Al. For all (2 ^ 0, w± € C'^{lRt x IRxx x 

A2. There is a constant C which does not depend on ^2 such that: 

\dxiW±{t,xi,^i)\ < C, y{t,xi,^i) GUJ+UUJ-. 


A3. There is a function which is continuous on a;+ and ui- such that w± tends to w± when 
^2 0 uniformly on the compact subsets of a;_|_ and a;_. 

A4. For all t* > 0 and x\ € IR, the limit 


lim w[ 


{t,xi,^i) (resp. liin w^{t, xi,^i)) 


exists and is denoted w!^{t*,xl) (resp. Wj.{t*,xl)). 


The characteristics related to (15.31) are the classical trajectories ipf {xj, ^i) = (t)) solution 

to: ^ 

f 1 a^i^(0)=X7 

1 ^ ^ ^ ^ ^ . (5.6) 

[A^f(t) = -a, e?(0)=6 

Due to the localization of the support of the initial data w/,+, only positive initial momenta > 0 
have to be considered in (15.61) . Therefore, the classical trajectories will reach the interface {.^1 = 0} 
at the time t* = ^. Moreover, due to the condition wj^- = 0, the behavior at the interface is 
described by the characteristics corresponding to the upper band only. To be more precise, by 
differentiating the map s 1 —>• w±(s,xj^(s),^(*"( 5 )), one gets by using (15.31) : 


W+(s,xj^(s),,^i+(s)) 


ds \ w-{s,xf{s),^f{s)) 


= -r(e+(s))M 


w+(s,xj^(s),,c^(s)) 

w_(s,x+(s),C+(s)) 


+ F{s) 


(5.7) 


where 

e^w = (e?(s),6), M = 


1 -1 

-1 1 


and F{s) = 


2^1+(s) 

|^+(s)| V ^xiW_(s,xj^(s),C+(s)) 


Applying the Duhamel formula to (15.7F it follows: 


w-{t,xt{t),^tit)) J V Ull,-{xi,^l) J Jo 

(5.8) 

Using (El), a direct computation shows that there is a constant C which does not depend on ^2 
such that: 


Vs, t e fR, \ f r(f±(^))d/r| <C. 

J S 


Moreover it holds 


lim / T{i^{s))ds = 

?2->0 Jo 


0 if 0 < t < U 
P At>t* 


(5.9) 


(5.10) 


In addition, since the lower band is not occupied initially and most of the band-to-band transfer 
occurs at the time t*, it holds: 


VO < t < t*, lim dxiW-(t, xj(t)) = 0 . 
^2—>-0 


(5.11) 


In order to simplify the presentation, the proof of (15.111) is postponed until the end of the present 
discussion. Consider first the case t < t* in (El- Using assumption A2 and equations (1^ and 
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(I5.11I1 . the dominated convergence theorem can be applied to the second term at the r.h.s. of 
equation (15.811 . which leads to: 


Jo 


(5.12) 


Combining assumption AZ with equations (15.1011 and (15.1211 , one can take the limit ^2 —> 0 in dSH) 
and obtain: 


f w\(t,xi +t,ii - at) \ ^ ^ wi^+{xi,ii) \ 
\ w^_{t,xi+t,^i - at) j \ Wi-{xi,^i) ) 


(5.13) 


To obtain (|5.13L we have used the limit below: 


lim (a:+(t),^+(t)) 

q2—>-U 


[xi + t, ^7 — at) if 0 < t < t* 

[xi + 2t* — t,^i — at) a t > t* 


which follows from the explicit formula of the classical trajectories solution to (15.611 . Then, by 
passing to the limit t —^ C in (j5.13p . it follows from assumption AA: 


( w'l{t*,xl) \ / wi,+ {xi,^i) \ 

wL{t*,xl) J Wl-{xi,^l) J 


(5.14) 


where xl = xi + t*. The case t > t* follows the same line with the difference that (I5.12|l has to be 
replaced by: 


Jo Jo Jt* 

The arguments used to deduce (15.1211 can be applied here to show that the first integral at the r.h.s 
of equation (15.1511 tends to 0 when ^2 —>■ 0. Moreover, using assumption A2 and equation (I5.9|l . the 
second integral at the r.h.s of equation (15.1511 is bounded by the quantity C{t — t*), which tends to 
0 when t ^ t* (here C is a constant which does not depend on ^ 2 )- Then, by taking successively 
the limit ^2 —0 and t —>■ t* in (15.8L we obtain: 


/ w\_{t*,xl) \ ^ -PM ( wi^+{xi,^i) \ 
w''_{t*,xl) ) wi-{xi,ii) ) 

Putting together (j5.14p and (I5.16L we get 

( w\{t*,x\) T \( w\{t\x\) \ 

\wUt\x\) ) \ T l-T )\wL{t\xl) ) 

where the transition probability T is given by: 

J. 1 -e-^^ 

2 


(5.16) 


(5.17) 


(5.18) 


and the constant (3 is defined in (15.511 . The system (15.1711 has the form of the solution of the well 
known Landau-Zener problem (see [3S]). Since xj and ^7 > 0 are arbitrary, equation (15.1711 is true 
for any t* > 0 and a;* S IR. 

We can now give the proof of (15.1111 . Differentiating (15.311 with respect to xi gives the following 
PDE for dx^W-: 


dtid, 


xiW-) - ^dx 


lei 


{dx^w-) - ad^^ (dx^w-) = T{^){dxjW+ - (dx^w-)) . 


(5.19) 


Consider {xi,^i) = [x^ (s),^;^ (s)), the lower band classical trajectory solution to (I5.6|l with an 
initial condition {xi,^j) such that (p^{xi,^i) = tpf {xi,^i). Then it holds 

dx,w-{t, x^ (t), {t)) = dx,w-{t, xt (t), (t)) 
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and it is enough to show that dx^W-{t,Xi (t),Ci (0) tends to 0 when ^2 tends to 0. Now, if one 
differentiates the map s (s)), one gets using (15.1911 : 

^ (s),5j"(s))) = r(C"(s)) w+(s, (s),(s)) - 52;iW-(s, a;]-(s), ^j"(s))) . 

By integrating the previous equation with the Duhamel formula, it follows: 

(9^,w_(t,x]-(t),Cr(t)) = f ^^^^'^'^T(C(s))dx,w+(s,x^(s),^^(s))ds (5.20) 

■Jo 

where we have used that u>/,_ = 0. Using assumption A2, together with the equations (15.91) and 
(15.1011 . we can conclude that VO < t < t*, the integral at the r.h.s. of (I5.20|l tends to 0 when 
6 ^ 0 . 

Remark 5.2. Although the system (15.1711 has the correct form, the formula (15.1811 for the transition 
probability is different from the Landau-Zener transition probability (13.311 obtained in section [3.21 
and which writes for our particular potential: 

T = e~T^\. (5.21) 


It will be verified numerically in section 16.31 that the band transmission corresponding to the 
effective model (15.311 is given by (I5.18|l whereas the band transmission corresponding to the model 
(15.211 is given by the correct transition probability (15.211) . This signifies that the approximation 
made in [5S] to obtain (15.31) from (15.21) is not correct. 


Remark 5.3. If the coefficient r is replaced by (iOD becomes a hyperbolic relaxation system 
(see [n, m) and the solution to (15.311 satishes w+ = W- when ^2 —>■ 0. Since / w+ + / W- = 1 
this leads to: 



(5.22) 


We will observe numerically in Figure IHl that (15.221) is true when the time is big enough. 


6 Numerical results 

In this section, the results provided by the surface hopping algorithm are compared with the 
reference level populations given by 63) where the solution u^{t) is computed numerically using 
an accurate method to solve the Dirac equation (12.1011 . In particular, it is verified numerically 
in section 16.11 that the spectral method is more accurate than the finite difference method. The 
image of the operator Il±(hD) is computed by using discrete Fourier transform (DFT) and Fourier 
multiplication. A comparison of the surface hopping algorithm with the models (15.21) and (15.311 is 
given in section 1^31 

6.1 Setup and quantum level simulations 

We suppose that the initial data Uj is such that its Fourier transform Uj satishes the relation: 

= f'\ox+m ( 6 . 1 ) 

where X+(0 is dehned in (12.911 and f^{^) is the Fourier transform of some function G L^{1R^). 
Such an initial data satishes the non-interference condition 63 - We remark that, like for the level 
population (14.11) . x+i^O can be replaced by X-i-(C) in (16.ip . 

If is bounded in L^{1R.^), using Lemma 2.3 in [TU] again, we obtain the following approxi¬ 
mation for the scalar Wigner transform of Uj: 

[u^}] = [f^] + 0(1) ( 6 . 2 ) 
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in V {rI X [Rl \ {^ = 0})) when /i ^ 0, where [/'^j = By plugging (H^]) in 

dUl), we obtain the following asymptotics for the initial value of the diagonal terms of the Wigner 
matrix: 

<+ = [f] + o(l), wl_ = o(l) (6.3) 

in V X (fR| \ {^ = 0})^ when h —>■ 0. In the present section, the initial upper level function 
is an h-scaled Gaussian wave packet: 


\x — x 


f^{x) = {'Kh) 2e 2h 


1^ , 

h 


with center Xg € R^, momentum ^g € R^ and norm = 1. Its h-scaled Fourier transform 

and its Wigner transform can be computed explicitly. Indeed, 




{-Kh) 


1 le-Col^ ,^n-( 

2 g 2h ^ h 


where Vm G L'^{R?), = h ^Tu{h Moreover, 




h 


In this section, the parameter h is equal to its physical value given by (1^ . The effective mass of 
the electron is given by 

m = 0.067TOe 

where me is the mass of the electron, the Fermi velocity is taken as vp = 10® m.s“^ and, having 
in mind the simulation of devices of size equal to hundreds of nanometers as in [ISJUl], we will 
take the reference length equal to L = 500 nm. Then, the numerical value of the parameter h 
corresponding to formula (12.711 is: 

/i = 3.4557 X 10"® 

which is small enough for the problem to be considered in the semiclassical regime. The simulation 
domain is equal to 

fl = [—IOa/Zi, lOVh] X [—5'\/h, bVh] ■ 

In the next section, we solve the Dirac equation (12.1011 for different choice of the potential U by 
using the time-splitting spectral method (TSSM) presented in Appendix |21 


6.1.1 The Klein tunneling 


The potential is equal to U = wol^i^o which corresponds to a Klein step. For this potential, we 
compare the TSSM and the Finite difference time domain method (FDTD) in [T^. The center 
and the momentum of the initial Gaussian wave packet are chosen as 

Xo = {-5Vh,0), fo = i(l,0). 

The height of the Klein step and the simulation time are respectively 

vq = 2|^o| and tf = IsVh. 

For the TSSM, the number of grid points are given by 


and for the FDTD by 


Ni = 126, N2 = 66 
Ni = 850, N2 = 426. 
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For the two methods, the time step size is 


At = 


All 

75 


where Axi and Ax 2 are the mesh sizes in the xi and X 2 directions respectively. In the present 
section and in the following section, the supscript h will be omitted for the solution u^{t,x) to 
(12.101) and for the initial condition u^{x). Using the discrete Fourier transform (DFT) (IA.5I) . resp. 
its inverse, to approximate the Fourier transform, resp. the inverse Fourier transform, we get the 
following approximation of the projectors Ii±{hD)u{t'^ ,Xj): 

= (6.4) 


where is the Fourier coefficient defined in (ESI and the discretization is the same as in 

Appendix E Then, using formula (14.111 . the approximation of the level populations 

of the Dirac equation is given by: 

^7±(i”) = l|n±(/ii^)«”ll2 (6.5) 

where \l±{hD)u'^ is defined by (16.41) and for u = {uj)j^j, uj € C^, we have 


\l = J2 ^2;iAx2 . 

jej 


( 6 . 6 ) 


Similarly, the initial condition ui{xj) defined by dSH is approximated by the formula: 

= Th --V E /"(efe)x+(/ia)e*«^"U J€J. 

J (&i-ai)(&2-02) ^ 


We remark that in the previous formula, we used the exact value of the Fourier transform f^{^k) 
instead of the DFT The initial data for the TSSM with the parameters described above is 

represented in Figure [TJ 

In Figure [H we depict the evolution with respect to the time t" of the level populations 
provided by the TSSM and by the FDTD. The curve with title Upper level, resp. 
Lower level, refers to the plus sign, resp. minus sign, in (16. 5|) and the curve with title Total refers 
to the total population Hu" 112- Initially, the charge is carried completely by the upper level, then 
non-adiabatic transfer occurs at time t = 0.28 and the charge is almost all transferred to the lower 
level. Moreover, we remark that the TSSM is more accurate than the FDTD. Indeed, the first 
method conserves the total mass (|A.8I) whereas for the second method the total mass decreases 
at the hopping time (it was shown in [T^ that the quantity conserved by the FDTD is not the 
mass but a related functional). To reduce this mass loss, the number of spatial points has to be 
chosen big enough which increases the CPU time of the method. In particular, the CPU times 
corresponding to the simulations of Figure[2]are the following: 493.5268 s for the FDTD, 1.3961 s 
for the TSSM. 

We remark in Figure[3]that the transition occurs when the charge reaches the Klein step which 
is a classically forbidden region for the Upper level. The band transition to the Lower level makes 
the Klein step an allowed region for the particles, which can tunnel in the region cci > 0 with a 
probability almost equal to 1. This is the Klein paradox, see |32) . 


6.1.2 Case U = axi, a > 0 


In this section, the potential is given by t/ = axi, a > 0. The Dirac equation (I2.10|) is solved using 
the TSSM. The center and the momentum of the initial Gaussian wave packet are taken to be 


4 = (-5v^,0), Co = (1,0). 
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Figure 1: Initial condition for the TSSM: |M/(a;)p (top) and, in a smaller region, Re ((M/) 2 (a;)) 
(bottom). 
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Figure 2: Time evolution of the level populations with a Klein step potential: comparison of the 
TSSM and the FDTD. 
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Figure 3; For different times and with respect to the space variable: representation of the position 
density of the projection of the wave function u solution to (I2.10|) where U = uola;i>o and the 
initial condition is given by (lO) . The upper half corresponds to |n+(/iZ))u"p and the lower half 
to |n_(/i£>)u"p where the projectors are computed using (16.4|) . The solution u is computed with 
the TSSM. 
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The level position density of the solution u{t,x) to (12.1011 is supported around {x^ (t), (t)), 

solution to: 

/ a;±( 0 ) = 4 

\ ^^±(t) = -a(l,0), C^(0)=eo 

The solution of the above problem can be computed explicitly. It is given by: 

x±(t) =x*T(|t-i*|,0), tit)=Co-ait,0) (6.7) 

where t* = is the time such that = 0 and the non-adiabatic transfer occurs. The point 

X* = x~^{t*) = Xq + (t*,0) is the point where the hopping occurs. The coefficient a is chosen such 
that the potential at the hopping point is equal to U{x*) = vq where vg = This leads to: 


^0 — (^o)i 

- 

(4)i 

The simulation stops at time: 

tf = ISa/Zi 

and the number of space points are given by 


Ni = 160, iV2 = 80. 


The time step size is equal to 


At = 


Aa;i 

TT' 


The initial data (uj)-, the Lower level and Upper level populations and the Total 

population are computed as explained in section Rl.l.ll The time evolution of the level populations 
is represented in Figured The numerical hopping is observed around the time t = 0.39 which is 
accurate enough compared to the predicted value t* = 0.3919 given by (EZD. Contrary to the case 
of the Klein step potential, a significant part of the charge stays on the upper level. 


The predicted value of the position corresponding to the hopping is x* = (9.7976 x 10“^, 0). 
We remark in Figure [S] that the transition occurs when the charge reaches x*. As for the Klein 
step potential, the band transition to the lower level allows the particles to tunnel in the region 
xi > xl- However, for the potential considered here, we can observe that the part remaining on 
the upper level is reflected. This could have been predicted from equation (EH). Indeed, for t > t*, 
the classical flow corresponding to the upper level (plus sign) moves to the left with respect to x* 
whereas the classical flow corresponding to the lower level (minus sign) moves to the right with 
respect to x*. 


6.2 The surface hopping algorithm 

In the present section, the potential U is equal to U = olx\ where a = 15. The initial condition is 
as described in section EH The center and the momentum of the initial Gaussian wave packet 
are taken as 

4 = (-5v^,0), Co = (1,0). (6.8) 

The diagonal terms w^{t,x,^) of the Wigner matrix defined by (14.3|) are computed using the 
surface hopping algorithm presented in section E) Then, for different values of the parameter h 
and of the simulation time tf, the surface hopping level populations ±(Z/) are computed from 
w\{tf ,by using (14.611 . The results provided by the surface hopping algorithm are compared to 
the level populations P^^^^{tf) of the Dirac equation computed as explained in section lH.l.ll 

For the level populations P^l^^^^tf), the Dirac equation (j2.10l) is solved using the TSSM in the 
simulation domain 

fl = [—IIa/Zi, llVh] X [— Sv^, bVh] ■ 
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Figure 4: Time evolution of the level populations for the potential U = axi, a > 0. 


Using equation (IA.8I) . we can choose a time step size which depends only on the simulation time 
f/ as follows: 

At = t//250. (6.9) 

The number of space points are given by: for h = 10“^ and h = 10“^ 


Ni = 250, N2 = 126, 


for h = 10-3 


iVi = 300, N2 = 150 , 


for h = 10"^ 


Ni = 850, N2 = 426. 


For the level populations ±(t/), the initial sampling is chosen as follows. We consider a uniform 
J X J discretization of the domain of the x variable: 


[ai, 6 i] X [ 02 , & 2 ] 

where oi = 02 = ccq — bVh and bi = 62 = Xq + By/h and a uniform K x K discretization of the ^ 
variable domain 

[Ci,di] X [C2,d2] 

where ci = C 2 = ^0 ~ and di = ^2 = Co + ^Vh. The grid points Xj, j = 1,.. ., and 
Cfc, k = 1 ,..., are ordered such that 
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t=0 t=0.1807 t=0.3976 t=0.5422 t=0.7591 



Figure 5: For different times and with respect to the space variable: representation of the position 
density of the projection of the wave function u solution to (I2.10I1 where U = axi and the initial 
condition is given by (EH). The upper half corresponds to |n+(h_D)u"p and the lower half to 
|n_(/ii7)u”p where the projectors are computed using (I6.4I1 . The solution u is computed with the 
TSSM. 
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Then, we determine the minimal numbers of points of the x variable and of the ^ variable 
such that 

^ > 1 - tou, ^ > 1 - tol5 

i=i fc=i 

where Aa; = , A^ = = ^2^£2. and tola,, tol^ are well chosen tolerances. The 

phase space points {xj,^k), j = 1, ■ • ■, fc = 1 ,..., are ordered such that: 

W’^ [f] (xi,a) > • • • > [f^] {xn^n,,^n^n,) 

and we determine the minimal integer N such that: 

N 

Y. w'^ [f^] {xt, ik)Ax^Ae > 1 - tol 

k=l 

where tol is a well chosen tolerance. This gives the first step of the algorithm of section S) Indeed, 
the initial sampling is the set: 

{{xk, Cfej +) G X fR| X {—, +} ; k = 1 ,..., N} 

where, using equation (1131), the associated weights Wk € IR can be approximated by: 

Wk = w'^ [f^] (xfe,a). 

In the present case, the hopping surface S defined by (13.21) is equal to: 

S = {{x,0&M*; 6=0} 

and the second order derivative of the function s i—>■ |'^^(s)P, defined by the characteristics solution 
to (12.121) . is equal to: 

(|e±(s)p)" = 2a2>0. 

Therefore, for the potential considered here, the points of extremal gap are all minimas. 
Moreover, the classical flow (x^(t), solution to (I2.12|) is given by 

x^it) = x±[ 6^(<) =^-a^(l>0). (6.10) 

In equation (16.101) . the formula of the momentum is explicit and the hopping transport step can 
be simplified. Indeed, for /c = 1,..., A^, if 0 < < 6’trajectory (xfe(<), ^k{t)), 0 < t < tj 

defined by (14. Sp will pass through a point (a:J,6) G R at the time t* = and non-adiabatic 
transfer occurs. For such a k, the weight is changed such that for t > t* 

Wk(t) = (1 - T*)wkit*) 

and a new particle is created on the lower band with index I > N. For the new particle, the 
associated weight is such that for t > t* 


wi{t) = T*Wk{t*). 

In the above equations, the transition rate T* is equal to 

T*=T{xUl) 

where T{x,^) is given by (15.211) . The surface hopping level populations defined by (14.61) are then 
given by: 

N M 

^iA.+(6) = ’ Psh,-{tf)= Y Mtf)Ax^A^^ . 

k^l /^AT+l 
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h 

PL.+{tf) 

Psl.+ {tf) 

PL.-{tf) 

Psh.-{tf) 

CPU dir (s) 

CPU sh (s) 

10 -^ 

7.48058 X 10-^ 

9.07411 X 10-"^ 

0.9251938 

0.9092579 

10.1206 

0.1400 

10 -^ 

8.98146 X 10-^ 

9.06980 X 10-^ 

0.9101854 

0.9093010 

10.1766 

0.1360 

10 -^ 

9.06998 X 10-^ 

9.06980 X lO-"" 

0.9093002 

0.9093010 

14.7849 

0.1320 

10-4 

9.06978 X 10-^ 

9.06980 X 10-^ 

0.9093015 

0.9093010 

198.3884 

0.1320 


Table 1: At time tf = 0.13 and for different values of the semiclassical parameter h: level popula¬ 
tions obtained by the Dirac solver and the surface hopping method. 


For all the tests, the size of the sampling grids are equal to: 

J = K =16. 


The sampling tolerances are taken as: 

tol = 10“®, tola; = tol^ = 10“^ X tol. 

For such parameters, the number of particles obtained numerically for the initial sampling is equal 
to iV = 6981. 

For the simulation time tf = 0.13 and for different values of the semiclassical parameter h, the 
level populations obtained by the two methods are listed in Table [TJ We remark that the reference 
lower level population -{tf) increases when h increases. This is due to the fact that for larger 
values of h, the transition process is slower and, at the time tf, the post-transition relaxation 
observed in Figure [3 has no yet happened. We notice that, for h smaller or equal to 10“^, the 
surface hopping level populations Psh±{tf) ^re almost constant with respect to h. This can be 
explained by the fact that, for the potential considered, the transition rate depends only on the 
variable which does not depend on h when considered on the sampling points. The column 

CPU dir, resp. CPU sh, denotes the CPU time required by the Dirac solver, resp. the surface 
hopping method, to complete the simulations. The time required to choose the initial sampling is 
not taken into account in CPU sh. It appears that the surface hopping method is very interesting. 
Indeed, when h decreases, the numerical cost of the Dirac solver increases whereas the surface 
hopping method provides very accurate results for an almost constant and much smaller CPU 
time. 

We depict in Figure [5] the absolute error of the level population corresponding to the upper 
level: 

£+{tf) ■■= \PL,+ itf) - Pk+{tf)\ (6-11) 

with respect to the semiclassical parameter h. We remark that the surface hopping level populations 
converge numerically when /i —>■ 0 to the level populations obtained at the quantum level. 

For h = 10“^, the time evolution of the level populations provided by the Dirac solver and 
the surface hopping algorithm are depicted in Figure [71 The curve with title Upper level dir, 
resp. Lower level dir, refers to resp. The same is true when dir is replaced 

by sh. The total population of the Dirac equation, curve titled Total dir, is obtained as explained 
in section 16.1.11 and the surface hopping total population, curve titled Total sh, corresponds to 
Psh + Psh.-{t)- We observe that the numerically obtained surface hopping total population 
is conserved. For a smaller CPU time, the time evolution of the level populations provided by 
the surface hopping algorithm agrees well with the one provided by the Dirac solver. Indeed, for 
the two methods, the charge is initially carried completely by the upper level, then non-adiabatic 
transfer occurs at time t = = II = 0.0667 (time required by the classical flow to reach the 

crossing set = 0} starting from the momentum of the initial Gaussian wave packet) and the 
great majority of the charge is transferred to the lower level. The CPU time required to complete 
the simulation is 14.7849 s for the Dirac solver and 5.6803 s for the surface hopping method. For 
the second method, the CPU time is bigger than the time indicated in Table [TJ This is due to 
the fact that the surface hopping curves in Figure [7| are obtained by repeating the surface hopping 
algorithm described in the present section for the sequence of times tf = 0.13n/100, I < n < 100. 
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Figure 6: At time tf = 0.13, logarithmic plot of the absolute error (16.111) of the level population 
corresponding to the upper level while varying the semiclassical parameter h = 10“^, p = 1,2,3,4. 
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Figure 7: For h = 10 ^ and U = axi, a = 15: time evolution of the level populations provided by 
the Dirac solver and the surface hopping algorithm. 
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6.3 Simulations using the models of [25] 


In this section we will compare numerically the models (j5.2p and (15.31) with the two dimensional 
version of our surface hopping algorithm. The comparison will be performed with h = 10“^, 
tf = 0.13 and with the same potential as is sectionThe initial condition is the 2D version of 
the initial condition in section EH In other words: 

w%^o = W^[f% w%=o = 0 


is replaced by 




... 1 («i-(a;n)i)^ (€i-(«n)i)^ 

{irhy^e - ^ -- 


w_(0,xi,^i) = 0. 


( 6 . 12 ) 


The functions w± and Wi depend on ^ 2 - However, since ^2 plays the role of a parameter, this 
dependence is not written on the l.h.s. of the equations in (16.121) . The center and the momentum 
of the initial Gaussian wave packet are given by (16.81) . 


6.3.1 The two dimensional surface hopping algorithm 


The a: 2 -independent solutions to the surface hopping algorithm presented in section 0] are obtained 
by replacing the system (I2.14|) with its two dimensional version: 

dtw± ± - ad^^w± =0, ^ 0. (6.13) 

More precisely, equation (16.131) is used for the time evolution of w± as long as the classical trajec¬ 
tories solution to (15.61) are away from the hopping surface {^1 = 0}. Then, the hopping 

transport is described as follows. 

Consider an initial set of sampling points 

{{xi,k, Ci,fe) +) S Mxi X X {—, -|-} ; fc = 1,..., N} 

with associated weights wt & M given by: 

Wk = w+{0,xi,k,ii,k) ■ 


For t > 0 small enough: 

(a:i,fc(t), Ci,fe(i)) = ‘ft ixi,k,ii.k), Wk{t) = Wk 

where, using (15.6|) . we have Ci,fc(0 = — Then, for fc = 1,..., A^, if 0 < < tf, the classical 

trajectory is such that ^i,fe(t*) = 0 at the time t* = and non-adiabatic transfer occurs. For 
such a k, the weight is changed such that for t > t* 


Wk{t) = (1 - T*)wk{t*) 


and a new particle is created on the lower band with index I > N. For the new particle, the 
associated weight is such that for t > t* 

wi{t) = T*Wk{t *). 


In the above equations, the transition rate T* is equal to: 

* 

T* = . 


The surface hopping level populations at the final time tf are then given by: 

N M 

wi{tf)^xiA^i (6.14) 

fc=l l=N+l 

where Aii and A^i are the mesh sizes in the xi and directions respectively. 
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6.3.2 Comparison with the models in |25| 

For equations (Ol) and (lOl) . the simulation domain is: 


(a:i,Ci) G 


— llVh, llv^ 


(Co)i 


atf — 5Vh, (^o)i + 5v^ 


They are solved with periodic boundary conditions on a uniform grid with 500 points in the xi- 
direction and 500 points in the ^i-direction. The time step size is chosen such that the condition 
of stability of the upwind method is satisfied [52]. More precisely, we take: 


At = 




Equations (Ol) and (15.31) are solved using a time-splitting method where the free equation (without 
source term) is solved using a dimensional splitting: the free problem is splitted in two one¬ 
dimensional problems and each one dimensional problem is solved using the one-dimensional second 
order upwind MUSCL scheme (see[22]). The source term is integrated in time using a RK2 method 
for equation (15.21) (in order to preserve the second order accuracy) and exactly for equation (15.31) . 

For equation (Ea, we take: 

Wi\t=o = 0 

in addition to the initial condition (16.121) . We remark that, since the MUSCL method is written 
for real valued functions, the third equation in (15.21) has to be splitted in two equations, one for 
the real part of Wi and one for its imaginary part. 

The level populations provided by the asymptotic model are defined by: 

Pam,±{t)= w±{t,xi,^i)dxidCi 

Jm? 

where w± is the solution to (15.21) and the level populations provided by the effective model are 
defined by: 

PL.,±(.t)= w±{t,xi,^i)dxid^i 

where w± is the solution to (15.31) . For the level populations ±{'tf) defined in (I6.14L the initial 
sampling is chosen as follows. We consider a uniform J x K discretization of the domain of the 
(a:i,5i) variables: 

[a, 6] X [c, d] 

where a = (a:Q)i — bVh, b = (xq)! -I- bVh and c = (^o)i — 5^/h, d = (^o)i + bVh. The grid points 
(a:i,j, ^i,fe), j = 1, ■. ■, J, k = 1,..., K are ordered such that: 


W+(0,Xi,i,^ip) > . . . > W+{0,XijxK,il,JxK) 


and we determine the minimal integer N such that: 

N 

iy+(0, ccl.fc, ^pfc)AxiA^i > 1 - tol 

fc=i 


where Axi = A^i = and tol is a well chosen tolerance. Then, the initial sampling is the 
set: 

{(a:i,fe, +) G X X {—, -|-} ; fc = 1,..., N} . 

For all the tests in this section, the number of grid points for the two dimensional surface hopping 
algorithm is: 

J=100, K=im 
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Figure 8: For h = 10“^, ^2 = 10“^ and U = axi, a = 15: time evolution of the level popula¬ 
tions provided by the asymptotic model (15.211 and the two dimensional surface hopping algorithm 
presented in section 16.3.1 1 


and the tolerance is 

tol = 10"® . 

We depict in Figure H] the time evolution of the level populations provided by the asymptotic 
model (15.21) and the two dimensional surface hopping algorithm. The curve with title Upper level 
am, resp. Lower level am, refers to Parn,+ {^)i resp. The total population, curve titled 

Total am, corresponds to P^m+i^) + The same is true when am is replaced by sh. 

The surface hopping curves are obtained by repeating the surface hopping algorithm described in 
section 16.3.11 for the sequence of times tf = 0.13n/500, 1 < n < 500. The behavior is the same as 
in section EH the charge is initially carried completely by the upper level, then hopping occurs at 
time t = 0.0667 and the great majority of the charge is transferred to the lower level. The time 
evolution of the level populations provided by the asymptotic model fits well the one provided by 
the two dimensional surface hopping: the surface hopping algorithm validates the model (15.21) . 

Figure El shows that non-adiabatic transfer is only partially recovered by model (15.31) . We 
depict the time evolution of the level populations provided by the asymptotic model (EH) and 
the two dimensional surface hopping algorithm where the transition probability T* is replaced by 
(15.181)(15.51) . The curve with title Upper level em, resp. Lower level em, refers to Pem,+ {'t)^ resp. 
Pem-{^)- The total population, curve titled Total em, corresponds to Pern,+ {^) + Pem,-{^)- The 
same is true when em is replaced by sh. As in FigureEl hopping occurs at time t = 0.0667. However, 
using the simplified model (15.31) . only about a half of the charge is transferred to the lower level 
which is close to the transition probability given by (15.181)(15.51) and substantially different from 
Figure El 

In Figure lTUl we verify numerically that the transition probability corresponding to the effective 
model (15.31) is given by the formula (15.181) . The curve with title Trans eff refers to Pem-itf), the 
population provided by the effective model EH on the lower level and at the final time, for 31 
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Figure 9: For h = 10“^, ^2 = 10“^ and U = axi, a = 15: time evolution of the level populations 
provided by the effective model (15.311 and the two dimensional surface hopping algorithm presented 
in section 15.3.11 where the transition probability T* is replaced by (I5.18|)(I5.5D . 
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Figure 10: For h = 10“^, ^2 = 10“^, tf = 0.13 and U = axi, a = 15: numerical verification of the 
formula (I5.18|) for the transition probability corresponding to the model (15.311 . 


different values of the constant /3 distributed uniformly on the interval [0, 5]. The constant [3 given 
by (15.51) is made arbitrary in (15.181) by replacing the coefficient r appearing in (15.31) by f = 

The curve with title Trans th is the representation of the coefficient T given by (15.181) for the 
same values of p. We remark that the two curves are very close which validates the limit ^2—^0 
performed in section 0 for the solution to (15.311 . 


Appendix A Time-splitting spectral method (TSSM) 


The solution u{t, x) to (12.611 is computed on the domain: 

Vt = [ai,6i] X [02,62] 

using a time-splitting spectral method as in [T3]. For r = 1,2, we choose the spatial mesh size 
lS.Xr = in the r direction for a given integer N^. Define a uniform grid 

= a-l-(jiAxi,j2Aa;2), j&J (A.l) 

where 

J = {j = e I 0 < ji < Ai , 0 < j 2 < N2} 

and a = ( 01 , 02 ). For a given time step size At, let u” denote the numerical approximation of 
u{t'^,Xj) at the time t” = nAt, n > 0. Then, is computed from u" by decomposing the 

problem (12.611 in the two sub-problems 

ihdtu = [—ihaidxi — iha2dx2] u (A. 2 ) 


and 


ihdtU = Uu. 


(A.3) 
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The free Dirac equation (IA.2I1 is solved using a spectral method in space and exact time integration, 
whereas equation (IA.3I1 can be integrated exactly on To discretize (IA.2I) . we introduce 

the trigonometric interpolant of w. 

(A.4) 


A^iA^2 


fcG/C 


where u[t,Xj) = u(t,Xj). In equation (IA.4I) . we have 




N\ N-2 


No 


K = {k = ih, fc2) e z" I - ^ < fci < ^< fc2 < ^} 


and 


= 27r 


ki 


bi — ai 62 — 0,2 


For a given function / : —>■ C^, fk is the discrete Fourier transform (DFT) of / defined by 


h = Yf^ 




(A.5) 


3^J 


where fj = f{xj). When / = is a sequence with fj € C^, the DFT of / is the sequence 

fk defined by equation (IA.5I1 . Inserting (IA.4I1 into (IA.2I1 . and using the orthogonality of the set 
|g*5fc-(a:-o)^ fc € /C} with respect to the scalar product of one gets the following system of 

ODE 

For any to € IR, the exact solution to the previous equation is given by: 

where 

M{S,0 = . 

Applying an inverse discrete Fourier transform, one obtains the following expression for the ap¬ 
proximation u(t)j of u(t,Xj): 

1 




keK 


We remark that using the eigenvectors (EH), the matrix B{f) can be diagonalized as follows: 

B{o = Pif)D{OPifr 

where 

D{f) = diag(|C|,-|CI) , Pif) = (x+(C),X-(0) 

and therefore 


M{S,0 = = 


cos(5|C|) -isin(5|^|)(^i - i6)/l^l 

-isin(5|C|)(^i -b*C2)/ICI cos((5|^|) 


. (A.6) 


The Strang splitting is the second order method constructed as follows: solve the first subproblem 
(IA.2I) over only a half time step of length ^. Then, we use the result as data for a full time step 
on the second subproblem (IA.3|) and finally take another half time step on (IA.2I) . This leads to 
the following method: 


N 1 N 2 




fcG/C 


u*/ = 


(A.7) 


= 


1 


A^iA^2 




fee/c 


32 





















Remark A.l. It follows directly from the unitarity of the DFT and of the matrix M{S,^) given 
by (IA.6I1 that the TSSM (IA.7I1 conserves the discrete total charge, i.e.: 

||u ”||2 = ||m °||2 , Vn>0 (A.8) 

where the norm ||.||2 is defined in (16.61) . Moreover, when the potential is equal to a constant 
U = Uo, the Dirac equation (12.61) admits the following plane wave solutions: 

which are integrated exactly by the TSSM (IA.7I) it k G IR? satisfies ^ = ^k' for some k' G tC. In 
the case of the Schrbdinger equation, an analysis of the stability of the TSSM was performed in 
[B] for initial conditions which are close to plane waves. 

Remark A.2. In the applications, the solution of the Dirac equation (12.61) moves away from 
the simulation box. Therefore, in addition to the periodic boundary conditions provided by the 
spectral method, we use absorbing boundary layers at the edge of the simulation box, see e.g. [5B]. 

Acknowledgements A.F. acknowledges Clotilde Fermanian Kammerer for discussion about sec¬ 
tion [XHl 
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